Nonlinear diffusion effects on biological population spatial patterns 



(N 

o 

(N 



Eduardo H. Colombo 1 and Celia Anteneodo 1 ' 2 
1 Department of Physics, PUC-Rio, Rio de Janeiro, Brazil 
2 National Institute of Science and Technology for Complex Systems, Rio de Janeiro, Brazil. 

Motivated by the observation that anomalous diffusion is a realistic feature in the dynamics of 
biological populations, we investigate its implications in a paradigmatic model for the evolution of 
a single species density u(x,t). The standard model includes growth and competition in a logistic 
expression, and spreading is modeled through normal diffusion. Moreover, the competition term 
is nonlocal, which has been shown to give rise to spatial patterns. We generalize the diffusion 
term through the nonlinear form dtu(x,t) — Dd xx u(x,t) v (with D,v> 0), encompassing the cases 
where the state-dependent diffusion coefficient either increases (y > 1) or decreases (y < 1) with 
the density, yielding subdiffusion or superdiffusion, respectively. By means of numerical simulations 
and analytical considerations, we display how that nonlinearity alters the phase diagram. The type 
of diffusion imposes critical values of the model parameters for the onset of patterns and strongly 
influences their shape, inducing fragmentation in the subdiffusive case. The detection of the main 
persistent mode allows analytical prediction of the critical thresholds. 
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I. INTRODUCTION 

Pattern formation in population dynamics has been 
studied both experimentally and theoretically. In exper- 
iments, the dynamics of insects and bacterial colonies, 
amongst others, have been observed From the the- 

oretical viewpoint, mean- field descriptions render macro- 
scopic or mesoscopic approximations to describe the be- 
havior of such complex systems. To take into account 
spatial inhomogeneities, one may construct a partial dif- 
ferential equation that rules the temporal evolution of 
the population density u(x,t), a function of the spatial 
position x and time t. Within this family of models, 
a standard one was first introduced by Fisher It 
consists of a reaction-diffusion equation, taking into ac- 
count growth and competition in the usual logistic form. 
Recently, a generalization of Fisher equation (FE) was 
introduced pHa], namely, 



d_ 

dt 



u(x,t) = DV 2 u(x,t) +u(x,t)(a-bJ[u(x,t)]), (1) 



where D, a, b are positive parameters and J is a func- 
tional of the density embodying nonlocality: 



J[u(x, t)] = / f(x,x')u{x',t)dx'. 
Jn 



(2) 



In the particular case f(x,x') = 8{x — x'), the original 
(local) FE is recovered. The introduction of the non- 
local form of the competition term is motivated by the 
consideration that products released in the environment 
by the individuals may either harm or support neigh- 
bors' growth. Interestingly, this nonlocal component was 
shown to give rise to the formation of steady spatial 
patterns [8|]. Diverse variants have been studied before. 
As influence functions f(x, x'), square and smooth forms 
[1, Q symmetric or not TC| have been considered. Non- 
locality in the reproduction rate [||, dimensionality 



and fluctuation effects [l2[ have been investigated too. 
In all those cases, however, spatial spread was described 
by normal diffusion. Meanwhile, there are indications 
that the spreading of biological populations is not due 
to purely random motion but influenced by the density, 
either to favor or to avoid crowding |13l - tl5| . Hence dis- 
persal is guided by a state-dependent diffusion coefficient 
rather than by a constant one. 

An important class of generalized diffusion equations is 
constituted by the porous media equation d t u = d xx u v , 
with v > [l6| . It is ruled by a state-dependent diffusion 
coefficient, proportional to u v ~ x , hence embracing the 
cases where the coefficient either grows or decreases with 
the density u. The generalization of Arrhenius law [l7j . 
the performance of thermal ratchets [lS] under this kind 
of diffusion have been studied. The nonlinearity leads to 
anomalous diffusion either superdiffusion for v < 1 
or subdiffusion for v > 1, recovering normal diffusion 
when v = 1. Microscopically, high density regions can 
slow down (y < 1) or intensify (y > 1) individual dis- 
placements, as a consequence of homophilic/phobic be- 
haviors that rule the dynamics of self-diffusion favoring 
or not the mobility among other individuals. While v < 1 
reflects a reaction to sparseness (with high diffusion coef- 
ficient where the density is low) , on the contrary, v > 1 is 
associated to immobilization in poorly populated regions. 

We will analyze the effects of nonlinear diffusion on 
pattern formation by considering the one dimensional 
generalized FE with nonlocal competition 
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u(x, t) = D——^u v (x, t) + u(x, t)(a — bJ[u(x, t)l) . (3) 
ox z 



Since alternative forms of the functional f(x,x') do not 
yielded substantially different results [8j, we will restrict 
our analysis to the case where f(x,x') is constant for 
x—w < x' < x+w, and zero otherwise, namely f(x, x') — 
-^Q(w— \ x — x'\), where is the Heaviside step function. 
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II. RESULTS 

Typical longtime patterns are shown in Fig. [TJ Numer- 
ical integration of Eq. ([3]) was performed by means of a 
standard forward-time centered-space scheme with inte- 
gration time step dt < 10~ 3 and width of spatial grid cells 
dx < 0.1. As initial condition we considered small am- 
plitude random perturbations either above the null state 
or around the nontrivial homogeneous solution. We set 
periodic boundary conditions, with periodic domain size 
L = 100. 




FIG. 1: Longtime patterns obtained from numerical integra- 
tion of Eq. d3j) , with a = b = 1, L = 100, D = 0.1, w = 10 
and different values of v indicated on the figure, in (a) linear 
and (b) logarithmic scales. In the inset, the longterm profiles 
in the full grid are shown. Fhe profiles are plotted for t = 200, 
but they remain unchanged after t ~ 100. 

Notice that, while the number of peaks is not affected 
by changing j/, the form of the patterns becomes sub- 
stantially different. By increasing v, the width (inverse 
concavity) of the crests increases and the density at the 
valleys decreases, such that for v > 1 disconnected re- 
gions can arise. 

Figure [2] shows the time evolution for v = 4, start- 
ing with small random values of the density u(x, 0). It 
rapidly increases for all x towards the level correspond- 
ing to the homogeneous solution, uq = a/b (t < 10), 
while patterns develop. After t = 100 no substantial 
changes are detected at the crests. Between successive 
crests, the density tends zero (exponentially fast with 
time). This fragmentation or clusterization process (2pj 
yields isolated population groups (clusters). Therefore 
fluxes between clusters are eliminated in the long time 
limit. This phenomenon is crucial in connection with 
"epidemic" spreads within a population. 




FIG. 2: Time evolution of the density profile obtained from 
numerical integration of Eq. © with a = b = 1, L = 100, 
D = 0.1, w = 10 and v = 4, represented for different times 
t indicated on the figure, in (a) linear and (b) logarithmic 
scales. The thick line corresponds to t — 200. 



A. Stability analysis 

To determine the stability conditions we follow the 
standard procedure of considering a first-order pertur- 
bation around the homogeneous solution u = a/b: 



u(x, t) = uq + e exp(ikx + Xkt) 



(4) 



where e is the perturbation initial amplitude, fc the wave 
number and is the exponential rate of temporal behav- 
ior. Substituting Eq. (jH) into Eq. (|3]) gives the dispersion 
relation 



-vDu^k 2 



sin(wfc) 
wk 



(5) 



that generalizes the one obtained by Fuentes et al. [9(. 
Defining the nondimensional rate Afc = Eq. ([5]) can 
be rewritten in a single parameter form as 



A fe = -p{wk) 2 - 



\{wk) 



wk 



with (3 = 



a wr 



(6) 



Negative Afc means relaxation back to the uniform state. 
Figure [3] depicts the dispersion relation in a typical case 
where A^ can take positive values allowing instability 
growth. 

Even if the analysis at short times does not guaran- 
tee the later evolution towards a stationary state, the 
mode with largest growth rate k* (absolute maximum 
of the dispersion relation A^. vs k) could play a crucial 
role. This mode will excite other wavelengths though the 
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FIG. 3: Dispersion relation = \k/a vs scaled k (solid line), 
for P — 5 x 10 -4 . The dotted line corresponds to the term 
— sin(wk) / (wk) and the dashed line to the zero line, drawn 
for comparison. In the inset, we show the position of the 
first maximum ko as a function of j3 = vDu^" 1 /(aw 2 ). The 
vertical dotted line indicates the instability threshold. 



coupling nonlinear term, however, if they are damped, 
k* will remain selected and its harmonics will shape 
the patterns. Substitution of the Fourier series expan- 
sion u(x,t) — J2kL-oo c k(t) exp(ifcx) into Eq. ((3]), when 
v — 1, leads to the following evolution equations for the 
Fourier coefficients c k [H 



dc k 
~~dt 



-Dk c k + ac k ~by j c m c k ^ m . (7) 



These equations are highly coupled through the last non- 
linear term. If v ^ 1, there will be still an additional 
nonlinearity in the first term of the righthand side, any- 
way let us consider the case where the first term is very 
small, allowing the existence of unstable modes. The am- 
plitude of the mode corresponding to the uniform state, 
Co, grows with rate a until stabilization, as observed in 
numerical simulations, e.g., in the example of Fig. [2] the 
level Mo = a/6 = 1 is attained at times of order 1/a. 
The mode with largest initial (positive) rate quickly de- 
velops and keeps dominating at intermediate timescales. 
Notice in Fig. [2] an almost perfect sinusoidal profile at 
time t ~ 20. If a unique mode contributes to the sum in 
Eq. ([7]) , it grows with rate given by Eq. ([5]) until stabiliza- 
tion while the remaining modes will be dumped. Actu- 
ally a set of undamped harmonics, characteristic of each 
value of v, also persists to shape the profiles. Typical 
Fourier spectra for the long-term patterns are shown in 
Fig.@] Although we do not have a rigorous mathematical 
proof, we will see that numerical results indicate that the 
dominant persistent mode, defining pattern wavelength, 
is the fastest growing one at short times. 

Afe possesses infinite local maxima located at k n , n — 
0,1,.... Since the absolute maximum is the first one, 
then k* = k (see Fig. [3]). In the inset of Fig. [31 k 
(numerically obtained) is plotted as a function of j3. For 
sufficiently small /3, A k is dominated by the last term in 



FIG. 4: (Color online) Fourier spectra for the longtime pat- 
terns shown in Fig. Q] The vertical dotted line indicates the 
position of the dominant mode ko- 



Eq. ©, yielding k = 9 /w with 9 ^ 1.43 tt. Fig. H 
shows that the dominant mode of long time patterns is 
in good accord with fco- 

Perturbations to the homogeneous solution vanish if 
Afc < 0. Since sin(wfc) is bounded, then the instability 
condition A k > implies 



(3 < (wkY 



(8) 



For the mode with largest growth, considering the ap- 
proximation k Q — 9q/w ~ 1.437r/u;, one has the instabil- 
ity condition 



P 



< 6»o" 3 ~ (1.43tt)- 



(9) 



This is equivalent to requiring the positivity of the first 
maximum. Notice in the inset of Fig. [3] that ko ~ 
1.437r/iu remains a good approximation in the whole in- 
stability region, below the threshold (vertical dotted line 
in the inset). Beyond this point the maximum becomes 
negative, hence the homogeneous solution recovers its 
stability for any wavelength. Equation © defines the 
critical value j3 c = 9q 3 for the onset of patterns. 

As intuitively expected, on the basis of the homoge- 
nizing role of diffusion, the inequality in Eq. indi- 
cates that the diffusion constant can not exceed a limiting 
value for the perturbation to depart from the homoge- 
neous state. In accord with Eq. (JSJ), in the limit D — > 0, 
patterns are also observed. Then, diffusion is not a nec- 
essary ingredient for the onset of patterns but has a role 
in pattern shaping. Figure [U[a) shows the density pro- 
files that emerge for different values of D in the normal 
case v = 1. For D = patterns are noisy due to the 
lack of the smoothing effect of diffusion and the ampli- 
tude is less uniform but the wavelength I is well denned. 
Moreover, between bumps, the density tends to zero as 
in the subdiffusive case of Fig. [5J The width of the bump 
at zero height, 2xo, is also well defined. Beyond fluctua- 
tions, the results are robust, at least under the two types 
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FIG. 5: (Color online) (a) Longtime patterns obtained from 
numerical integration of Eq. ([3]) up to time t = 100, with 
a — b — 1, v = 1, w = 5, L — 100 and different values of 
D indicated on the figure. In the inset, the profiles in the 
full grid are shown. In (b) the time evolution for D = 10 _J is 
shown in logarithmic scale. Times are indicated on the figure. 



of initial conditions analyzed. In the absence of diffu- 
sion, the steady state must verify u(x)[a — bJ{x)] = 0, 
then either u{x) vanishes or its integral within the inter- 
val (x — w, x + w) must adopt the constant value a/b. The 
former case requires that the null solution becomes sta- 
ble in some regions. The latter requires that each cluster 
does not see the neighbors and any point in the cluster 
must be influenced by the full cluster. This means that 
2xq < w < £ — 2xo, which is verified in numerical exper- 
iments. 

For D — 10~ 5 , patterns are still noisy at t = 100, 
but are expected to smooth out at much longer times. 
In this case also the density between crests goes to zero 
exponentially fast with time, as can be seen in Fig.(5jb). 
In the case of the figure, the clusterization occurs for 
D < 10~ 3 while for larger values of D not only the crests 
but also the valleys stabilize in a finite value as time goes 
by. Then clusterization occurs for D below a threshold 
value only. It is noteworthy the resemblance of the noisy 
profiles with those observed in experiments with bacteria 
[4J. But here clusters arise naturally, without imposing 
absorbing nor zero flux boundaries. 

If D 7^ 0, Eq. © predicts the existence of a minimal 
value of the interaction range w required for pattern for- 
mation, with all other parameters kept fixed. This crit- 
ical value does depend on the kind of diffusion through 
the factors v and Uq _1 . Notice also that for v ^ 1 there 
is influence of uq too, which is absent in the normal case 

(!/ = !). 



According to the hypothesis that ko is the characteris- 
tic wavenumber of the stationary pattern, and taking into 
account that boundary conditions are periodic (i.e., an 
integer number of wavelengths must accommodate into 
the size of the system L), the number of maxima m is 
given (in average) by 



knL On L L 
m = — = — - ~ 0.715- . 

2f 2fw w 



(10) 



Even in the cases when Eq. (I10|) gives an integer value, 
it is expected to furnish the number of peaks observed in 
the average. In practice, depending on the initial condi- 
tions, the crests come out and grow accommodating its 
number approximately to the rounded value of m. Fluc- 
tuations in the effective m are larger the larger m or the 
further is m from an integer value. For instance in the 
example of Fig. [TJ instead of seven, eight crests are ob- 
served in some realizations, being m ~ 7.15. 

These observations can be verified through numerical 
integration of the evolution equation. In Fig. [SJ we show 
the number of maxima m of the patterns as a function of 
w together with the theoretical prediction given by Eq. 
(ITU1) . An excellent agreement between theoretical and 
numerical outcomes can be observed. Then, in good ap- 
proximation, the dominant wavelength only depends on 
the relation L/w independently of the remaining parame- 
ters. However, these parameters participate in condition- 
ing pattern upraise through the critical value of /3, given 
by Eq. ©, and may also influence pattern shape. In par- 
ticular, in the average, m does not depend on u, as can 
be observed in the example of Fig. [TJ The wavenumber is 
preserved in the limit D — > 0, even if in approaching this 
limit more and more modes become unstable (i.e., more 
local maxima of A^ become positive) , but the global max- 
imum remains approximately the same. Then, diffusion 
is not necessary for pattern formation, nor influences the 
characteristic wavelength. 




w 



FIG. 6: Number of maxima m as a function of the nonlo- 
cality width w, for a = b = 1, L = 100 and D = 0.01, and 
different values of v indicated on the figure. The solid line 
corresponds to the approximate theoretical relation Eq. (|10|) 
and the symbols to the outcomes of numerical simulations. 
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monotonic behavior of the critical curve is broken when 
uq < 1. Thus, for low values of w there is also an upper 
critical value of v for the onset of patterns, but for large 
w, patterns occur for any v. 



B. Patterns shape 

Although the characteristic mode does not depend on 
u, its amplitude does. This is shown in Fig. [9] where the 
amplitude, Au = u max — Umin obtained from numerical 
simulations is represented as a function of w for different 
values of v. 



FIG. 7: Number of maxima m as a function of v for two 
different values of w indicated on the figure, with a = b = 1, 
L — 100 and D = 0.01. Dashed lines correspond to the 
theoretical prediction given by Eq. (|10[1 . with the additional 
condition ([9]), defining a critical value v c beyond which no 
patterns occur (we set m = in such case). 

Notice also that condition © indicates that, although 
approximately the same number of maxima is expected 
independently of v, there are critical thresholds v c be- 
yond which no patterns occur. This is illustrated in Fig. 
[7] for the case a = b = 1, for which v c — w 2 /(D6q). 

Differently from normal diffusion, uo = a/b is also de- 
terminant of pattern formation. How the phase diagram 
is altered by changes in uo is illustrated in Fig. HI The 
shadowed area represents the region where no patters 
emerge when uq = 1. For other values of uq, only the 
frontier is shown. For uq > 1, the critical curve increases 
monotonically with v, such that only small v (superdif- 
fusion) allows pattern formation for a given interaction 
range w and remaining parameters fixed. However, the 
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FIG. 8: Phase diagram of pattern formation in the vw plane 
, for wo = 1, following Eq. (|9}. Shadowed is the region where 
no patterns arise. The lines show the frontier of the phase 
diagram for other values of uq = a/b. Patterns emerge for 
in > i» c (above the critical line). Parameter values are L = 
100, D = 0.01 and a = 1. 
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FIG. 9: Pattern amplitude Au = Umax — u m in as a function 
of the interaction width w, for a — b = 1, L — 100, D = 0.01, 
and different values of v indicated on the figure. The dotted 
lines are a guide to the eyes, the vertical ones indicate the 
critical value predicted by Eq. |9j. 

For vanishing w we recover the local case f(x,x') = 
S(x — x') in which no patterns emerge, as supported by 
numerical simulations. In agreement with Eq. @, there 
is a critical value w c , at which the amplitude vanishes. 
Notice the abrupt decay of the amplitude at the criti- 
cal value. This threshold was not detected in previous 
works dealing with normal diffusion possibly because of 
the range of parameters used. For instance in Ref. [8j, 
w c /L would be of the order of 10~ 3 . The critical value 
w c decreases with v, indicating that a shorter influence 
range w is required when the dispersion passes from sub- 
diffusive to supcrdiffusive. Then superdiffusion favors 
pattern formation and the amplitude of the patterns is 
larger. For w = L/2 (or its multiples), the nonlocal 
term becomes J[u(x, t)] — J(t), that follows the equa- 
tion dJ/dt = (a — bJ)J. Then, in the long time limit 
J — > uq, implying that the homogeneous state should 
also be attained in this extreme case. 

Furthermore, v can have strong effects on patterns 
shape. As depicted in Fig. [T] subdiffusion (v > 1) in- 
duces fragmentation (clusterization) . Solutions that van- 
ish outside a finite support are typical of nonlinear subd- 
iffusion [Tj|. In the opposite case v < 1 (super-diffusion), 
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the effects are not so striking concerning pattern shape, 
for moderate values of the diffusion coefficient. The re- 
gion between crests assumes larger values the smaller 
v. Fragmentation also emerges for any kind of diffusion 
when D is small enough (as discussed in connection to 
Fig. [5]) or also if w becomes large enough (not shown). 
The shape of the clusters depends on v. Their amplitude 
decays and their width increases as v increases. 

It is noteworthy that the distance between crests 
(wavelength), i — L/m ~ lAw, is larger than the interac- 
tion range w, however if the cluster size 2xo(w) is large 
enough, there can be influence of one cluster over the 
two neighboring ones. When clusters are disconnected, 
2^0 ;$ i — w — 0.4k; means that one cluster does not 
influence the neighbors. Otherwise they do, even if dis- 
connected. 



III. FINAL REMARKS 

Nonlinear diffusion is expected in the spreading of bi- 
ological populations rather than normal diffusion, hence 
motivating the introduction of a state dependent diffu- 
sion coefficient, as in Eq. Q. We have shown how pattern 
formation is altered in the presence of anomalous diffu- 
sion. Moreover, in all cases, the initially fastest growing 
mode remains selected at longer times. This observation 
allows to obtain theoretical predictions that we verified 
through numerical integration of the evolution equation. 



Then, it is clear that diffusion is not a necessary in- 
gredient for the onset of patterns, nor has impact on the 
characteristic wavelength, that depends only on the in- 
teraction range w. Furthermore, diffusion imposes a crit- 
ical threshold of the model parameters for pattern forma- 
tion. The type of diffusion regime has impact on patterns 
shape, even if the characteristic mode is kept unchanged. 
An important qualitative change in the shape of patterns 
occurs mainly for v > 1, in which case fragmentation of 
the population is induced. This effect is also observed for 
very small diffusion constant D and/or large interaction 
width w. The occurrence of fragmentation may have im- 
portant consequences in diseases dissemination and other 
spreading processes triggered by contacts between indi- 
viduals. Superdiffusion (v < 1) facilitates pattern forma- 
tion, which can occur even for shorter interaction width 
w and manifesting larger amplitudes than in normal dif- 
fusion. 

Beyond the initial motivation of introducing nonlinear 
diffusion to the nonlocal FE, we uncovered aspects that 
apply also to the normal diffusion case previously studied. 
The identification of the main mode selected at long times 
allows to perform analytical predictions, which may be 
extended to tackle other variants of the model. 
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